Pattern formation in a predator-prey system characterized by a spatial scale of 

interaction 
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We describe pattern formation in ecological systems using a version of the classical Lotka-Volterra 
model characterized by a spatial scale which controls the predator-prey interaction range. Analytical 
and simulational results show that patterns can emerge in some regions of the parameters space 
where the instability is driven by the range of the interaction. The individual-based implementation 
captures realistic ecological features. In fact, spatial structures emerge in an erratic oscillatory 
regime which can contemplate predators' extinction. 
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The study of spatial aspects of population dynamics, 
with the analysis of patterns and the mechanisms deter- 
mining spatial correlations, started with a seminal work 
by Moran [1]. Afterwards, theoretical studies highlighted 
how different synchronizing mechanisms might interact 
to produce spatial patterns [2]. At present, a central is- 
sue in population ecology is whether the predator-prey 
interactions can be considered among these mechanisms. 
In fact, spatial correlations between preys and predators 
can be observed in real populations and a clear empirical 
evidence for this phenomenon can be found, for example, 
in a system involving predatory beetles and larval flies as 
their preys [3] . Generally, these systems have been theo- 
retically described using Lolka-Volterra inspired models 
which are capable of generating diffusion-driven instabil- 
ities. The analysis of these reaction-diffusion models ar- 
gued that development of spatial patterns is possible only 
under some general conditions [4] regarding the values of 
the diffusion coefficients (predator disperses faster than 
the prey) and regarding the type of the growth functions 
and predator functional response [5] . 

In this paper we are interested in showing how differ- 
ent and more realistic mechanisms can generate spatial 
inhomogeneities which are not directly produced by dif- 
fusion phenomenon. Standard Lotka-Volterra equations 
describe a two-species system where both species can co- 
exist with the population densities regularly oscillating in 
time. This result is similar to empirical observations and 
it is a remarkable prediction considering the simple math- 
ematical elegance which characterizes the model. With 
the aim of preserving the lightness and economy of this 
approach, we will introduce a simple modification of the 
Lotka and Volterra's original idea in a spatial version of 
their model. The new ingredient is a probability of in- 
teraction (predation) which becomes a function of the 
distance between individuals. An intuitive justification 



rests upon considering that the probability that a con- 
sumer meets a prey should be dependent on the relative 
distance between them [6]. The introduction of a spa- 
tial scale of interaction has been widely applied to model 
competition describing species coevolution in community 
ecology [6] . More recently it was used in studies of evolu- 
tionary theory which propose models for the emergence 
of polymorphism or sympatric speciation [7]. 

We consider a model characterized by a couple of equa- 
tions, one for the prey N(x, t) and one for the predator 
P(x,t). They describe diffusion in real space and the 
strength of the interaction in the nonlinear term is a func- 
tion of individuals' proximity [8]. We do not introduce 
a single interaction scale L because we consider effective 
ranges of interaction (the real meeting area can have a 
different relevance for the growth of predators and for the 
death of preys): 
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Predators consume the preys with an intrinsic rate a 
and reproduce with rate /?; r is the preys' growth rate 
and predators are assumed to spontaneously die with 
rate m. and Dp are the diffusion coefficients of 

preys and predators, respectively. This system presents 
two stationary and spatially homogeneous solutions: an 
absorbing phase N(x,t) = P(x,t) — and a survival 
phase N(x,t) = ^j^', P(x,t) — 2a r Li . With the aim 
of investigating the existence of solutions with spatial 
structure we make a stability analysis around N(x, t) 
and P(x,t) by considering small harmonic perturba- 
tions [4,8]: N(x,t) = N + A N exp[Xt + ikx}; P(x,t) = 
P + Apexp[Xt + ikx}. Their introduction into eq. 1 leads 
to a linear system which exhibits solutions when the de- 
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terminant equals zero. This condition results in the fol- 
lowing dispersion relation: 



k 2 / fc 4 

A(fe) = - — {D N + D P ) ± J — (£>j\ 



Dp) 2 



sin{kLi )sin(kL2) 



k 2 L ± L 2 



(2) 



This relation shows a symmetry in the interchange of 
the diffusion constants or the interaction lengths. Spatial 
patterns can emerge if the condition Re[X(k)] > is satis- 
fied. For this reason, L\ ^ Li is a necessary prerequisite. 
For the case Dn = Dp = 0, the condition obviously 
allows patterns formation. This fact proves that the in- 
stability is driven by the range of the interaction and is 
independent of the diffusion process. It is interesting to 
note that in this case N(t) = Y^ij=i a,j(t)5(x — jAzjv) 
and P{t) = Y^j=i bj{t)${ x — j^-xp) are solutions of 
the system [9], with ai(t) and bi(t) solution of rescaled 
classical Lotka-Volterra equations. From these results 
arise that our instabilities are not diffusion-driven, where 
Dn =/= Dp is a necessary, albeit not always sufficient con- 
dition for generating spatial patterns. It follows that it is 
natural to simplify our analysis taking Djy = Dp = D. 
Moreover, for L 2 = 2Li and by introducing the rescaled 

variables K = kL\ and A = \-fj-, eq. 2 reduces to: 



X(K) = -K 2 + 



DK 



■\f— sin 2 K cos K . 



(3) 



This relation is shown in Fig. 1. The onset of the in- 
stability can be identified by the values of the parame- 
ters for which the maximum of the curve becomes zero 
(X(K m ) = 0). We can compute it numerically and, for 
the original variables, we obtain: 



k r 



1.82759 



™ l > 12.5232. 



D 



(4) 



It is easy to extend this analysis to a two dimensional 
space. There the new relations read: \k\ m w 2.19535/Li 
and s/rmL\/D > 22.4228. 




FIG. 1. Dispersion relation X(K) for different values of the 
parameters (/i = yJrmL{/ D). 



Now we propose a microscopic discrete stochastic for- 
mulation of the model which allows us to describe the role 
of demographic fluctuations. This implementation intro- 
duces two relevant differences. The first one is due to the 
role of intrinsic stochasticity which causes internal noise. 
The second one is specifically related to the discrete na- 
ture of individuals which can generate threshold effects 
not present in a continuum description where every small 
amount of the density of population is acceptable, an as- 
sumption of continuity which is often unrealistic (atto-fox 
problem [10]). Individual-based lattice models for simi- 
lar Lotka-Volterra system have been studied in details in 
previous works [11,12]. In general, the introduction of 
stochasticity generates a system considerably richer and 
perhaps even more realistic. 

For the sake of simplicity our individual based model 
is implemented in a 1-dimensional system where the fol- 
lowing algorithm was carried out. Simulations start with 
an initial population of Pq predators and N preys, ran- 
domly located along a ring (periodic boundary condi- 
tions) of length equal to 1. The different processes of 
diffusion, reproduction and death are implemented se- 
quentially by randomly selecting an individual of each 
population (predator or prey). The selected action is re- 
peated for a number of times equal to the size of the corre- 
sponding population. When all the processes are carried 
out, a time step ends and the algorithm restarts. In de- 
tail, we are considering five processes: 1) diffusion, where 
a predator (prey) is randomly selected and moves some 
distance, in a random direction, chosen from a Gaussian 
distribution of standard deviation a. 2) predator repro- 
duction, with rate /3-/Vf; 2 , where JV£ is the number of 
preys which are at a shorter distance than L 2 from the 
predator at position x. 3) predator death, with probabil- 
ity m. 4) prey reproduction, with probability r. 5) prey 
death, with rate aP]^ , where P|J is the number of preda- 
tors which are at a shorter distance than L\ from the 
prey at position y. All the newborns maintain the same 
location as the parents. We evaluate N£ 2 and PJJ using 
periodic boundary conditions. If, in eq. 1, we measure 
time in units of the simulation time step, the coefficient 
D is related to the discrete model through D = a 2 / 2. 
Birth and death probabilities are the same in the con- 
tinuous and in the discrete model. A rigorous derivation 
which would show that the continuum field equations ap- 
proximating the discrete model correspond to the eq. 1, 
can be obtained by using Fock space techniques [8,12]. 

The agent based simulations produced a rich collec- 
tion of data which allow to explore the temporal and 
spacial behavior of our system. Irregular oscillations, 
which swing in a rather erratic fashion around an av- 
erage value, characterize the temporal evolution (sec 
Fig. 2). Other stochastic models display a similar behav- 
ior [11,12]. There, inevitable fluctuations tend to push 
the system away from the trivial survival phase and in- 
duce irregular population oscillations that almost resem- 
ble the deterministic cycles of the classical Lotka-Volterra 
model. These features well approximate the temporal 
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evolution of our simulations characterized by homoge- 
neous spatial distributions. For example, the behavior of 
the dominant Fourier component displayed by the time 
evolution of these data depend only on the parameters 
r and m, in a fashion consonant with a classical Lotka- 
Volterra system. Moreover, the presence of oscillations 
is independent of the population size, and therefore they 
persist in the thermodynamic limit [12]. The more the 
spatial solution is marked by clustering (lower L\ and 
D values), the more the irregularities of the temporal 
oscillations increase. Finally, we must remember that 
an important outcome of the introduction of intrinsic 
stochasticity is the possibility of predators' extinction. 
For obviously, when the number of predators becomes 
very low, a chance fluctuation may lead the system into 
a state with P(x, t) = 0. Therefore, asymptotically as 
t — > oo this state will be reached. 




FIG. 2. Temporal evolution. Top, populations with ho- 
mogeneous spatial distributions (Li = L2 = 0.1). Bottom, 
populations with modulated spatial distributions (Li = 0.05, 
Z/2 = 0.1). Others parameters are: a = 0.004, r — m — 0.5, 
a = /3 = 6.25 x 10~ 5 , P = N = 40000. 

Now we turn to the analysis of the spatial aspects of the 
distributions of the populations. In accordance with the 
results of the continuous model no regular patterns can 
emerge for interactions characterized by the same range. 
For L\ 7^ L2, clear spatial patterns are generated (Fig. 3). 
For an interaction equal to the ring dimension, the simu- 
lation is characterized by a spatially homogeneous occu- 
pancy. Decreasing the L\ value one peak appears, with 
the population concentrated in one region of the ring. 
For lower values of the interaction length, some regions 
of the ring are occupied forcing all the remaining areas, 
up to some range, to be nearly empty. This state corre- 
sponds to a sequence of isolated colonies (spikes) [8,13]. 
We also observe that high population density can gen- 
erate distributions where the colonies merge up, without 
loosing the ordered character of the modulation. We have 
deeply explored the case with L\ = 2L 2 . Even if 7^ L 2 , 
the number of peaks is generally equal for predators and 
preys and the tuning of the parameter L\ allows modu- 



lations of arbitrary wavelengths. Finally, for extremely 
short-ranged interaction, a noisy spatially homogeneous 
distribution appears (see Fig. 4). The reported periodic 
spatial patterns are stationary, with configurations char- 
acterized by a noise which increases in the neighborhood 
of the transition towards the homogeneous distribution. 
These outcomes obtained from simulations on the seg- 
ment (0, 1) can be extended into a 2-dimensional space 
where fluctuating clusters arranged on an hexagonal lat- 
tice can emerge. 

We can observe that for low D values and low popula- 
tion density disordered spikes can appear, independently 
on the values of L\, probably generated by the asymme- 
try between birth and death processes [14]. In fact, birth 
events only occur adjacent to a living organism, whereas 
deaths occur anywhere. This fact introduces a source of 
spatial correlation which can result in a reproductively 
driven cluster mechanism. Hence, under reproductive 
fluctuations and in the presence of a weak diffusion, indi- 
viduals can organize into clusters. These inhomogeneous 
configurations can be related with the patterns observed 
in discrete lattice simulations, where spatial structures 
can also be generated by traveling waves [12,15]. 
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FIG. 3. Spatial distributions. Top, homogeneous spatial 
distributions ( L\ — L2 = 0.1). Bottom, modulated distribu- 
tions (Li = 0.05, L2 = 0.1). We can note preytaxis. Same 
parameters as in Fig. 2, t=600. 

In the previous paragraph we show how inhomoge- 
neous spatial distributions can appear, depending on 
the parameters value. Now, we will try to charac- 
terize the transition towards these states (segregation 
transition). A proper order parameter is provided by 

N(t) 2 



q M = max q>0 



^2 exp[i27rq • xj(t)] , where the sum is 
3=1 . 

performed over all individuals j with their positions de- 
termined by Xj at a given time r. The transition from a 
homogeneous to an inhomogeneous distribution matches 
the jump of qm to an integer value, corresponding to 
the number of periodic clusters present in the space. In 
fact, if the space is homogeneously occupied qM — 1-4 
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and the segregation transition is characterized by the 
passage of qm from 1.4 to an integer value as soon as a 
modulation becomes dominant [13,8]. In Fig. 4 we show 
qm as a function of D and L\. For any value of the 
range of the interaction, a critical value of the diffusion 
coefficient (D c ) exists above which no spatial structures 
emerge. For r = m = 0.5 and L\ = 0.1 the analytical 
prediction gives D c — 4 x 10~ 4 , in good accordance with 
the discrete model where the first set of simulations for 
which qm = 1.4 presents a diffusion coefficient equal to 
D c 3.7 x 10~ 4 . Moreover, a critical value of L\ exists 
(L c ) for which the segregation transition takes place. For 
L\ > L c clusters appear and the value of L c obtained 
by the simulations is in accordance with the analytical 
result coming from eq. 4. Finally, another correspon- 
dence between the predictions of the continuous model 
and the Monte Carlo simulations exists. In fact, the con- 
tinuous description can even reproduce quantitatively 
the period of the patterns in the modulated distribu- 
tions. For example, considering the case displayed in 
Fig. 3, for Li = 0.05, the first relation in (4) tell us that 
the fastest growing mode is k = 36.6. This is in good 
agreement with 6 x 2ir = 37.7, which is the wavenum- 
ber of the configuration generated by our simulation. In 
fact, this wavenumber is the first immediately above the 
fastest growing mode of the mean-field description which 
is compatible with the periodic boundary conditions. For 
obvious, from the same equation in (4), we can obtain 
the general analytical relation for the number of peaks n: 
n w 0.29X7^ , which is compared with the Monte Carlo 
data in Fig. 4. 



In conclusion, the introduction of a finite-range in- 
teraction in a spatial Lotka-Volterra model, allows the 
description of a rich spatio-temporal dynamics charac- 
terized by regular spatial structures. This type of spatial 
behavior is a central issue in population ecology as, 
effectively, it is possible to record spatial correlations 
between preys and predators in nature. Our investiga- 
tion was carried out both analytically, using a suitable 
continuous model (mean- field approach), as well with a 
discrete model, by employing Monte Carlo simulations. 
The individual-based model shows the existence of spa- 
tial structures in an erratic oscillatory regime which can 
contemplate predators' extinction. The mean-field de- 
scription captures the essential features of the discrete 
model. We record quantitative correspondences for the 
period of the spatial patterns in the modulated distribu- 
tions, the value of the critical diffusion and the critical 
interaction length. From an ecological perspective, the 
emergence of spatial patterns in an erratic oscillatory 
regime are realistic elements generally absent from con- 
ventional approaches and might shed further light on 
issues of particular ecological relevance. 

We thank the Brazilian agencies CNPq for partial fi- 
nancial support and the school CSSS 2008, of the Santa 
Fe Institute, where this work was conceived. We are 
grateful to B. Perthamc for helpful discussions. 



3.5 | i i i | i i i i | i i i i | i i i i | i i i i | i i | i 
3 • • • • _ ! 



1.5 - 



0.005 0.01 0.015 0.02 0.025 0.03 
a 




0.01 0.026 0.064 0.16 0.4 



'-I 

FIG. 4. Top, qui as a function of a [L\ — 0.1, L2 = 0.2, 
r = m = 0.5, Po — No = 20000). Data are averaged over 50 
realizations. The dashed line stands for the analytical result 
a c ~ 0.028. Bottom, log-log plot. qM as a function of L\ 
(L 2 = 2L 1 ,<t = 0.004, r = m = 0.5, P = N a = 40000). Data 
are averaged over 20 realizations. The continuous line rep- 
resents the analytical prediction: n w 0.29/Li. The dashed 
line stands for the analytical result L c ~ 0.014. The structure 
functions are averaged over 20 time steps, starting at t = 580. 
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